arXiv:1509.05521v2 [math.NA] 3 Mar 2017 


NOVEL RESULTS FOR THE ANISOTROPIC SPARSE GRID QUADRATURE 


A.-L. HAJI-ALI, H. HARBRECHT, M. D. PETERS, AND M. SIEBENMORGEN 


Abstract. This article is dedicated to the anisotropic sparse grid quadrature for functions 
which are analytically extendable into an anisotropic tensor product domain. Taking into ac¬ 
count this anisotropy, we end up with a dimension independent error versus cost estimate of 
the proposed quadrature. In addition, we provide a novel and sharp estimate for the cardinality 
of the underlying anisotropic index set. To validate the theoretical findings, we present several 
examples ranging from simple quadrature problems to diffusion problems on random domains. 
These examples demonstrate the remarkable convergence behaviour of the anisotropic sparse 
grid quadrature in applications. 


1. Introduction 

This article is dedicated to the construction of anisotropic sparse grid quadrature methods, 
where we emphasize on the Gauss-Legendre quadrature. Anisotropic sparse grid quadrature meth¬ 
ods can be seen as a generalization of sparse Smolyak type quadratures, cf. [21], since they are 
explicitly tailored to the anisotropic behaviour of the underlying integrand. Taking into account 
these anisotropies leads to a remarkable improvement in the cost of the sparse grid quadrature. 

Usually, a sparse grid quadrature is described by some sparse index set and a sequence of 
univariate quadrature rules. For the sequence of univariate quadrature rules, we employ Gauss- 
Legendre quadratures with linearly increasing numbers of quadrature points. The index set is a 
priority chosen with respect to a certain weight vector, which incorporates the anisotropy, and a 
predefined approximation level. It is also possible to adaptively select those indices which provide 
the main contribution to the integral, see m- Such adaptive methods have successfully been 
applied in the context of random diffusion problems, see e.g. (T] |22l |26], in order to compute 
best iV-term approximations of the corresponding solution. However, the adaptive construction 
of index sets is computational expensive and only heuristic error indecators are available. Hence, 
it can not be guaranteed that the adaptively selected index set is optimal. Instead of choosing 
Gauss-Legendre points, a sequence of nested quadrature rules such as Clenshaw-Curtis or Leja 
type quadratures could be considered as well. While the number of quadrature points needs 
to be doubled for Glenshaw-Curtis quadratures in order to guarantee their nestedness, only one 
additional quadrature point is added for each consecutive member of the Leja sequence. Hence, 
based on the Leja sequences, a sparse grid quadrature can be constructed where only one additional 
function evaluation is required for each new multi-index in the sparse index set, cf. m- However, in 
contrast to Gaussian quadratures, the quadrature weights of the Leja sequence are not necessarily 
positive which yields that the stability constant of the quadrature might not be uniformly bounded. 
Moreover, Gaussian quadrature rules provide a much higher degree of polynomial exactness than 
Leja quadratures with the same number of quadrature points, which is particularly advantageous 
for smooth integrands. 

The main task in estimating the quadrature’s cost is the estimation of the number of multi¬ 
indices which are contained in the sparse index set. For the isotropic variant, the number of 
indices can easily be determined by combinatorial arguments, see e.g. uniES]. Things get more 
involved if one considers anisotropic, i.e. weighted, sparse index sets. In this case, to the best of 
our knowledge, only very rough estimates on the cardinality of the index set are known, although 
several estimates can be found in the literature, see e.g. [Ij. In fact, this problem is equivalent to 
the estimation of the number of integer solutions to linear Diophantine inequalities, see m and the 


2000 Mathematics Subject Classification. 65D30, 65C30, 60H25. 


1 



2 


A.-L. HAJI-ALI, H. HARBRECHT, M. D. PETERS, AND M. SIEBENMORGEN 


references therein, which is a problem in number theory, or to the calculation of the integer points 
in a convex polyhedron. Current estimates are not sharp and do not provide improved results for 
the cost of the anisotropic sparse grid quadrature in comparison with the anisotropic full tensor 
product quadrature. In this article, we prove a novel formula to estimate the cardinality of the 
sparse index set in the weighted case. This formula is much sharper than the other established 
estimates. 

A very popular application that requires efficient high-dimensional quadrature rules are para¬ 
metric partial differential equations. They are obtained, for example, from partial differential 
equations with random data by truncating the series expansions of the underlying random fields 
and parametrizing them with respect to the random fields’ distribution. As representatives for 
such problems, we shall consider here elliptic diffusion problems with random coefficients as well 
as diffusion problems on random domains as specific examples to quantify the performance of 
the anisotropic sparse grid quadrature. The resulting quadrature approach is very similar to the 
anisotropic sparse grid collocation method based on Gaussian collocation points which has been 
introduced in [231 El] • The collocation method interpolates the random solution in certain col¬ 
location points and represents it in the parameter space with the aid of polynomials. Thus, it 
belongs to the class of non-intrusive methods, cf. [1]. Instead of representing the random solution 
itself, the anisotropic sparse grid quadrature can be employed to directly compute the solutions 
statistics, i.e. its moments, and functionals of the solution. 

The remainder of this article is organized as follows. Section [^specifies the quadrature problem 
under consideration and provides the corresponding framework. The subsequent Section |3| is 
dedicated to the construction of the anisotropic sparse grid quadrature method. In particular, the 
main ingredients, i.e. the index set and the sequence of univariate quadrature rules are introduced. 
In Section 131 we provide corresponding error estimates with respect to the level of the anisotropic 
sparse grid quadrature and with respect to the cardinality of the index set based on the one 
dimensional error estimate for the Gauss-Legendre quadrature. Section [5l deals with the cost of 
the anisotropic sparse grid quadrature. In particular, we state here a novel estimate on the number 
of indices in the weighted sparse tensor product and provide a proof of this estimate. In Section|6l 
we consider three different applications: A pure high dimensional quadrature problem, a diffusion 
problem with random coefficient and a diffusion problem on a random domain. Finally, we state 
concluding remarks in Section [71 


2. Problem setting 

In what follows, let T := [—1,1]. The cr-algebra of Borel sets on T shall be denoted by B 
and v '■= dy/2 is the normalized Lebesgue measure on B. We define the product probability 
space (T°°,B°°,fi) of all sequence^ ip: N* — r, where xp = {V'rajn- Herein, B°° is the cr-algebra 
generated by the cylindrical sets and ^ is the corresponding product measure, i.e. /r(Ai x • • • x 
^rr.xrx...) = n:=iK^o for all m S N and Ai,..., Am & B. For an integrable function 
/: F°° —>■ K., we are interested in the efficient approximation of the integral 

(1) [ /d^- 

Jr°° 

The approach we shall present here is based on the construction of efficient quadrature formulas 
for a certain surrogate fm '■ T™ —>■ R of /. In order to define this surrogate, we make the following 
assumption. 

Assumption 2 . 1 . Let = E(r,T„) := {2 G C : dist(z,F) < t„}. We assume that f is 
analytically extendable into S(x) := Sn for an isotone sequence Tn —t 00, which measures 

the anisotropy of the function f with respect to the different dimensions. 

Now, given an anchor point xp G F°“, we define the surrogate fm of / as the projection of / 
onto the first m variables anchored at xp, i.e. 

fm{yi^ • ■ • ) 2/m) •= fivij • ■ • ) 2/m) V’m-l-l) V'm-l-2) • ■ 

^We make the convention N := {0,1, 2,...} and N* := {1, 2,...}. 
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The projection /„ is well-defined since / is analytic. Moreover, we know from Assumption 12.11 
that fm is analytically extendable into 

m 

■= X 

n—1 


and it holds that 
( 2 ) 


J /dli-y /m(y)2 ™dy 


< ce(m), c > 0 


with a null sequence e{m). 

In the sequel, we shall approximate 


(3) 


lfm- = 


^I^”^)/™:=y^/m(y)2-™dy 


by the anisotropic sparse tensor product quadrature. Herein, we set 

(l^”Vm)(y*) := y /m(2/n,y*)^, 

where y* := [yi,..., yn-i,yn+i, ■ ■ ■, 2/m] with the notational convenience y = (j/„, y*). It is evident 
that the precision of the applied qnadrature has to increase when m increases. Moreover, the cost 
usually scales exponentially with respect to m, which is referred to as the “curse of dimensionality”. 
Therefore, we have to keep track of the impact of the dimension m on the error estimates. 

We start by considering the univariate quadrature with respect to y^ n = 1,... ^m. To that 
end, we introduce the Wpoint Gauss-Legendre quadrature operator with quadrature points 
€ r and weights Wk according to 


N 

(Q^"Vm)(yn) 


k=l 


This quadrature operator is exact for all polynomials p S n 2 Ar-i := span{I, x,..., More¬ 

over, : C'(r) —>■ R is continuous with continuity constant 1. Hence, the quadrature error can 
be bounded by 


(4) 


|(i'”>-Q(”))/m(y*)|< inf 

pen2N-i 

< 2 inf 
pen2«' 


(|l^"n/m(y^) -P)| + - P)\) 

\\fmiyn)-p\\ 


lc(r)- 


Thus, in accordance with e.g. [H Chapter 7.8] and [T], the analytic extendability of fm guarantees 
that the A^-point Gauss-Legendre quadrature satisfies the one-dimensional error estimate 

(5) - Q^”^)/m(y;^)| < c(K„)exp(-log(«;„)(2iV- I)) ||/m(yn)||c(E„) 


with 

Kn = Tn + \/ 1 + and c(«;„) 
Here and in the sequel, we set 


4 


Kn- 1 


ll/m(y];)||c(E„) := max 1/^(0, y*)|. 


Note that it holds by definition 

4 

> 1 -I- r„, Kn > 2Tn and, thus, c(k„) < —. 

Tn 
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3. Anisotropic sparse grid quadrature 


We shall now introduce anisotropic sparse grid quadrature formulas which extend the original 
idea of Smolyak’s construction from |29) . To that end, we start by considering an increasing 
sequence of univariate quadratures Q^- with points and weights 


( 6 ) 


j = 0,1,2,..., 


where Ni < N2 < ■ ■ ■ and Nj —)> 00 for j — >■ 00 . 

Since we deal with multidimensional quadrature rules, we furthermore define tensor products 
of univariate quadrature rules by 

/ ra s 

Qn ® ® E • ■ • E { 

ki = l km=l ^i=l ^ 

Following the notation of [25] , we introduce for j G N the difference quadrature operator 

(7) Aj := Qj - Qj_i, where Q_i := 0. 

With the telescoping sum Qj = isotropic m-fold tensor product quadrature, which 

uses Nj quadrature points in each direction, can be written by 


(8) (g) • • • (g) AW (g) • • • (g) , 

l|a||oo<j 

where the superscript index indicates the particular dimension. Since the tensor product quadra¬ 
ture is convergent, we observe that 

(9) I(/,^)= ^ AW®...®AW/,^. 


Jm) 


The cost of applying a quadrature formula is measured by the number of quadrature points. In 
case of the isotropic tensor product quadrature operator ([8]) this number is determined by {Nj)^. 
Thus, this isotropic tensor product quadrature suffers extremely from the curse of dimensionality. 
The classical sparse grid quadrature, cf. [SI HU], can overcome this problem up to a certain extent. 
It is based on linear combinations of tensor product quadrature formulas of relatively small size. 
To define the sparse quadrature, we introduce as in [2l|23| for each approximation level q the sets 
of multi-indices 

X{q,m) := |o < a G N™ : ||a||i < 

and 

F(q,TO) := |o < a G ■ q — rn < \\a\\i < q'^. 

The sparse grid quadrature operator, cf. [a HU] HU], is then given by 

(10) A{q, m) := ^ A« 0 • • • ® AW. 

OL^X {q,m) 

An equivalent expression is obtained by the combination technique, cf. El, 

(11) A{q,m)= V (-l)«-ll“lb(' where Q„ := Q« ® • • • 0 Qi™’- 

A visualization of the set of indices X{q,m) is given in Figure (Tj 

The number of quadrature points used in or m is considerably reduced compared to the 
full tensor product quadrature. However, the sparse quadrature operator does not take into ac¬ 
count the fact that the different parameter dimensions are of different importance to the integrand 
fm- Indeed, the cardinality of the set X{q, m) is given by 

which still depends heavily on m. Thus, we assign a weight to each particular dimension and use 
a weighted version of the Smolyak quadrature operator. 
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Figure 1. The 21 indices contained in the sparse grid X{5,2) on the left and 
the 56 indices contained in X(5,3) on the right. 


Let w € denote a weight vector for the different parameter dimensions. We assume in the 
following that the weight vector is sorted in ascending order, i.e. wi < W 2 < ■ ■ ■ < Wm- Otherwise, 
we would rearrange the particular dimensions accordingly. We modify the sparse grid index sets 
X{q,m) and Y{q,m) in the following way, see also 


( 12 ) 

and 


X^{q, m) := < 0 < a e N™ : ^ < q 


n—1 


(13) 


Y^{q,m) 


0 < a e N” 


m 

: q - ||w||i < ^ anWn 

n—1 



With this notation at hand, the anisotropic sparse grid quadrature operator of level 9 € N is 
defined by 


(14) A^{q,m):= AW®...®AW 


which can equivalently be expressed as, cf. [M] . 


(15) A^{q,m)= ^ Cw(Q;)Qa, with c^{a) := ^ (-I)II^IIl 

aeY„(q,m) 

The formula (1151) can be regarded as the anisotropic combination technique quadrature. For the 
evaluation of this formula, we only need to determine the coefficients Cw(q:) and to apply tensor 
product quadrature operators of relatively small size. Thus, in order to compute the approximation 
to © with the anisotropic sparse grid quadrature (na, it is sufficient to evaluate the integrand 
fm on the anisotropic sparse grid 

’Avji.dl IR) .- X * * * X 

cx.G.Yw{q,m) 

Note that the sparse grid quadrature operator coincides with the anisotropic sparse grid 
quadrature operator (fT4l) for the special weight vector w = 1 := [ 1 , 1 ,...,!]. 

In Figured the indices of the weighted sparse grid A(i_ 2 . 5 )( 5 , 2 ) and of the weighted sparse 
grid 2 , 3 )( 5 ,3) are visualized. We observe that the number of indices is drastically reduced in 
comparison to the according isotropic sparse grids visualized in Figure [TJ 
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Figure 2. The 10 indices contained in the weighted sparse grid -^( 1 , 2 . 5 ) (5) 2) on 
the left and the 16 indices contained in -^( 1 . 2 , 3 ) (5, 3) on the right. 


The computation of the anisotropic sparse grid quadrature operator (fldll depends on the choice 
of the weight vector w and the sequence in ([ 6 ]). In view of the one-dimensional error 

estimate and to keep the number of quadrature points low, the sequence is chosen 

slowly increasing in accordance with 


(16) 




ia+ 2 ) 


Then, we can find an upper bound for the contribution of the difference quadrature operator 
Aj = Q^- — Qj_i for all j > 1 and for all functions /i: T —^ R which are analytically extendable 
into Si according to 

|A,/i|< |I/i-Q,-/i| + |I/i-Q,_i/i| 

<c(iti)(l + e-i°8(«^))e-'°s('‘^)^||/i||c(E0 
< (!l±Ic(ici)e-'°s(«^)^||/i||c(E,) < 2c(Ri)e-i°«(-^)^||/i||c(E,)- 

Kl 

For j = 0, the difference quadrature operator coincides with the function evaluation at the mid¬ 
point z = 0 of F which implies that 

(18) lAo/il = IQo/il = |/i(0)| < 

Analogously, it follows from ([SJ and (ITBl) that 


(19) |I/i - Q,/i| < c(ici)e-'°*5('^^)(^+i)||/i||c(E,)- 

Next, let us consider the multivariate integrand fm '■ F™ —>■ R which can be analytically extended 
into the region Due to the analyticity in a tensor product domain, it follows that the 

contribution of the tensor product of the operators Aj is bounded by the product of the one¬ 
dimensional contributions. Indeed, we obtain that 


(20) 


< sup |(ld®Ai 2 ) ® • • • 0 AW)/„(z) 


- ( n (2c(«:n))’ 


“in(l.“-) ]g-E™=ilog(« 


n—1 


'’W/niWcCEm) 
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with ||/m||c(Sm) |/m(z)|- In addition, we take the minimum in in order to ensure 

that the constant is 1 if = 0 in accordance with (1181) . 


4. Error estimation for the anisotropic sparse grid quadrature 

For the estimation of the quadrature error of the anisotropic sparse grid quadrature, we employ 
the following lemma. 

Lemma 4.1. Let {ipn}n G be a summable sequence of positive real numbers. Then, there 

exists for each <5 > 0 a constant c(5) independent of q > 1 such that 


( 21 ) 


+ 1) < c{S) exp{Sq). 


n —1 


Proof. Let 0 < (5i, ^2 < ^ be arbitrary such that ^ 1+^2 = 6 . From the summability of {ipn}n, it 
follows that there exists a jo = jo(^i) € N such that 


( 22 ) 


ifn< dl. 

n=io + l 

We now split the left-hand side in m into 

00 JQ 00 

(23) n (9'0n + 1) — (O'V'n + 1) {<11pn + !)• 

n—1 n—1 n—jo-\-l 

Then, the second factor can simply be estimated by 


00 z 00 \ 

{qifn -b 1) = exp f ^ log(gV'n + 1)) < exp((5ig). 

= '>.-. 4-1 ^ n—jo-\-l ' 


‘n.=jo + l 


The number of factors jo in the hrst product on the right-hand side of (1^ is fixed and depends 
only on the choice of di and on the decay properties of {ipk}k- Since jo is a fixed natural number, 
there exists for all <52 > 0 a constant c{Si,S 2 ) such that 


30 


Yliqipn + 1) < c((5i,d2)exp((52g)- 


n —1 


Hence, we obtain that 


{qifn + 1) < c((5i, ^ 2 ) exp(dg). 


n —1 

Since 0 < i5i, ^2 < <5 can be chosen arbitrary with the only limitation that -I- ^2 = d, the choice 
c{5) = infij+ij^i 0 ( 61 , 62 ) yields the desired estimate. □ 

With the above preliminaries and further assumptions on the summability of the sequence 
{7Vi}n, we are able to establish error estimates for the anisotropic sparse grid quadrature. 

Assumption 4.2. The sequence which describes the regions of analytic extendability of 

the function f, fulfills 

Tn P Cn^ 

for some r > 1 and a constant c > 0. Hence, the sequences {T~^}n and {(k„ — 1)“^}„ are 
summable. 

For the error estimation, we apply an identity which was used in |24) to bound the error of 
a collocation approach. Additionally, we exploit Assumption 14.21 to obtain an estimate which is 
exponentially decreasing in the sparse grid level q and does not depend on the dimensionality m 
at all. Therefore, we further notice that the integration operator I: C(Sm) —?> R is obviously 
continuous with continuity constant 1. 
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Lemma 4.3. Let the sequence of quadrature points be chosen as in (HU) and let the weight vector 
w be given by Wn = log(K„). Then, there exists for each (5 > 0 a constant c{S) independent of m 
such that the error of the anisotropic sparse grid quadrature is bounded by 

(24) \{1-A^{q,m))fm\ < c{6,T)e~'^^^~^^\\fm\\c{-E,n) 
with 

c{S,t) = 4:C{S)\\{T-^}n\\ii 

and c{S) denotes the constant from Lemma \4.1\ with respect to the summable sequence 
Note that the constant c{S, t ) tends to infinity as S tends to 0. 

Proof. In the same way as in |24j . the error of the sparse grid quadrature is rewritten, with the 
notation I = by 

m m 

(25) l-A^{q,Tn) = Y,R{q,n) (g) 

n=l k—n-\-l 


Herein, the quantity R{q, n) is defined for n = 1 by 


R{qA) 

and for n > 2 by 


R{q,n) 


n—1 

E 0<’® 



Q 


L(9-e 


n-1 
k = l 


Oik'^k 



Due to R{q, 1) = , we deduce for the first summand in (1^ that 


m \ 

k=2 ' 

The summands in (l?5|) with n > 2 can be estimated with (IT^l) . (1^01) and with the continuity of the 
integration operator by 


R{q,n) 


t('=) 


k—n+1 


< 




n ( 2 c(Kfc))' 


a_l ^ k—1 


min(l.afc) \ Y^Zl “fc log(Kfc) 


< c(k„) 


g-log(Kn)([(9-Efc;) akWk) lwr,\+l) -YAZl afclog(Kfe) 


' n—1 


n ( 2 c(Kfc)) 


min(l,afc) 


k^l 


ll/m||c(5 


With the choice Wk = log(Kfc) for all fc = 1,..., m, it follows that 
(R{q,n) (g) I«V„ 

^ k^n +1 ^ 


< c(k„) ^ 


g-g-I 2 k=i Q^kWk 




OiGXw2^..^_i(g,n — l) 


' k^l 
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It remains to estimate 


E 


n i^Kk)) 


min(l,afc; 




- E (n (2cK)) 

citGX^{q,m) ^ k—1 




The expression inside the product always equals 2c{Kk) except for the case ak = 0. Hence, it 
follows that 


E ( n ( 2 c(Kfc)) 

v.GX^{q,m) ^ k—1 


min(l,Q:fc) 


L—I 1^1 

L u;i J ^wm J 

\ min(ai ,1) 


< ( 2 c(ki))‘ 


E! ( 2 c(/Cm)) 


min(aTn,l) 


ai =0 


Ckrn —0 


< 


n 


/ 2 c{Kk) q 

V Wk 


+ 1 ) < c(S) exp(Sq). 


The last inequality holds since {2c(Kk)/wk}k is summable and, thus. Lemma l4Tl is applicable. 
Combining our findings, we obtain that 

m 

\{1- A^{q,m))fm\ < 11/^110(5:™) 




which yields the estimate (PTI) . 


□ 


Lemma 14.31 implies that the anisotropic sparse grid quadrature converges exponentially with 
respect to the level q. The convergence in Lemma 14.31 is nearly as good as the convergence of the 


anisotropic tensor product quadrature on level q, with [—h i] quadrature points in the n-th 
dimension. 

In addition, we provide an estimate on the quadrature error in terms of the number of multi¬ 
indices contained in the anisotropic sparse index set. Note that this is very similar to the analysis 
in From ®, we deduce that the error of the anisotropic sparse grid quadrature can be 

written as 


(26) 


I (l (o'! ^))/m I — 


q:€N^\Xw(<7,7Ti) 


With estimate (l20)) . we obtain 


\{l-A^{q,m))fm\ < E ( n ( 2 c(Kn)) 


min(l.a„) \ log(K„)o 


‘ll/m||c(S„) 


ct£N'^\Xw{q,m) ^ n—l 

O'G N'^ \ Xw ( O': in) 


with 


(k) := sup ( 2 c(k„))' 


[lin(l,Q;n) 


OcGN' 


n—1 


Due to the summability properties of {t„}„, the constant c{k) can obviously be bounded indepen¬ 
dent of the dimension m. It remains to estimate the sum in the above estimate which has been 
extensively studied in m- The following result from m is particularly useful for the considered 
situation. 

Theorem 4.4. Let the sequence w = {wn}n of positive, real numbers be ordered, i.e. Wi < wj for 
i < j ■ If there exists a real number j3 > 1 such that 


A^(w,/3) := E ■ 


1 


n—1 




(27) 


< oo. 
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it holds that 
(28) 


E 

olGN^\Xw{q,c 






We apply Theorem 14.41 with respect to the sequence = log(A„). Then, it follows from 
Assumption 14.21 that the conditions of Theorem 14.41 are fulfilled with /3 < r since then 


(29) 


E 


1 

glog(K„)/,S _ 2^ 


^E 


1 

(1 + cn’')i//3 - 1 


< oo. 


Of course, Theorem 14.41 is still valid when considering dimensions m < oo. In this case, even 
subexponential convergence rates can be proven which, however, depend on the dimension to, see 
m for the details. Especially, the appearing constants depend then on the dimensionality as 
well. Since we are interested in dimensionalities which might grow with the desired accuracy, it 
is reasonable to rely on estimates which do not depend on to. To that end, we conclude from 
Theorem 14.41 for all /3 < r that 


(w,/3) 

(30) 1(1-Aw (g, to)) /mI < c(k) - - - < c{w, 


with (w,^) := E^=i —. Due to (|29l) . the latter constant is bounded independently of 

TO. 


5. Cost of the anisotropic sparse grid quadrature 


5.1. A preliminary estimate on the cost. In order to find an error estimate in terms of 
the number of quadrature points, we additionally have to estimate the cost of the sparse grid 
quadrature method for a given level q. 

In the following, we establish a bound on the number of quadrature points used in the combi¬ 
nation technique formula (I15|l . This number is obviously bounded by 


m 

cost (Aw (g, to)) < ^ 

ct£Ywiq,rn) 


which might be a rough upper bound since some of the coefficients in (1151) may vanish and since 
some of the quadrature points usually appear repeatedly in m- Since Y^{q,m) C X^{q,m), 
cf. (fT^ and ([13]), we can further estimate 


(31) 


cost (Aw (g, to)) < ^ J([ 

CK.GY^{q,m) n—1 


-(a„ -I- 2) 


< ^ JJ(«n + l) 

cx.^Yvjiq,m) n—1 


< n(«"+i)^ 

cx.£Xw{q-,ni) n=l 


max I I (an H- 1) 

x^X^(q,m) 

n—1 




Note that there holds f3 € X^{q,m) for all 0 < /3 < a whenever a e X^{q^m). This property is 
often referred to as downward closedness of the index set. Hence, it follows that 

m 

#Aw(g, to) > max > 1= max 1 f (ci„-I-1). 

0<P<oc n—1 

As a consequence, the number of quadrature points in (IT^ with a sequence of univariate quadrature 
sequence, where the number of points are given by (HU), is bounded by 

(32) cost (Aw(g,TO)) < #Aw(g,TO)^. 

We remark that a similar bound for downward closed index sets has also been used in |9] in the 
context of a sparse adaptive collocation approximation. 

As indicated before, the estimate (1^ is not sharp. In fact, we see that cost (Aw(g, to)) depends 
rather linearly than quadratically on the number of multi-indices ^X-w(q,m), cf. Figures |4]-|6| from 
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the numerical examples. This is due to the fact that an exact representation for the cost is given 

by 

m 

(33) cost = E nc.. 

cx.GX^{q,m) n—l 

where denotes the number of quadrature points which belong to but not to for 
any i < a„. In our setting of the Gauss-Legendre quadrature, where the number of points is 
determined by (HID, this sequence is given by 

{Cn}n = {1,2,0,2,0,4,0,4,0,6,0,6,...}. 

Therefore, each summand in (1331) vanishes whenever any is an even number greater than zero. 


5.2. A sharp estimate on the anisotropic sparse index set. In order to complete the con¬ 
vergence analysis, it remains to estimate the number of indices in the set X-w{q,rn). To that end, 
we require the following lemma. 

Lemma 5.1. For L G N, m G N and S G M+, there holds the inequality 

L~1 m ^ m 

j—O n—l n—0 

with equality when d = 0. 


Proof. We prove the assertion by induction on L. For L = 1, we verify 

m +1 


n(»^+<5)= 


1 


1 


m-\-t + 5 


n ^ (n -b d -f 1). 


n—l n—l n—0 

Let the assertion be fulfilled for L. Then, we conclude for L + 1 that 

L 

■ Tl) (L T 5 + 77.) 


j=0n—l n—l 

' L -G m -\-1 + 5 


n—l 


777. + 1 


(L + 77 + 5) 


n—l 


1 


TO + 1 


m+1 


(L + n -b d) 


n—l 


1 


TO -b 1 


tt (-b + 1 + H -b l5). 


n—0 


□ 


The next lemma gives us a novel bound on the number of indices in X^{q^ to). 


Lemma 5.2. The cardinality of the set X^{q, to) in (|T^ . where the weight vector w = + 1 ,..., Wm] 
is ascendingly ordered, i.e. wi < W 2 < ■ ■ ■ < Wm, is bounded by 


(34) 


my N 

#Aw((?, m) < -b 1). 


Proof. The prove is performed by induction on m. For 777 = 1, the assertion is obviously fulfilled, 
since 



#Aw(<z,l)= ^ 1 

ai=0 



-bi. 
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Let US assume that (IM)l is true for m — 1. For m G N, the cardinality of X^{q, m) can be calculated 
by 

^ ^ (9 jWrmTn !)■ 

i=i 

Inserting the induction hypothesis yields that 


(35) 


j—0 n—1 

/ m —1 


nWr} 


n 1 

fc=l ^ 

71 — 1 / 

n 1 


—I 


+ 


^ n—1 
' m — 1 


nb 


n—1 


nwn 


1- Wrn -I m — 1 1 I q—jWrr 

E n - 

j—0 n—1 

E n(i- 

n—1 

/ L^^l m-l 

1+E n(i- 

\ j=l n=l 


JWr, 


JW-n 


nWn + q 


Focusing on the last term and since Wm 7 Wn foi' all 0 < n < m, we conclude that 


E n(i 

j=l n=l 


JWr, 


^E n(i 

i=l n=l 
m—1 

' n 


JWr, 


n 

n—1 


Wr, 


nWm + q 
-1 m-l 

E n 

i=l n=l 


n - 


Applying the previous lemma with L = and S = X - L leads to 

m— 1 


L m—1 L — 1 m — 1 

E n ^=E n “ n ^ 

j—l n—1 j—0 n—1 


Thus, we obtain that 


Lll^Jm-l 

E n(i 


JWm 


1=1 "=1 

Inserting this into (1551) finishes the proof. 


nWr, 


< 


L + 5 


m 


n—0 


mWr) 


□ 


Remark 5.3. (1) We would like to point out that estimate (1341) is sharp in the isotropic case, 

that is, for the weight w = 1. Moreover, the ordering of the weight vector is crucial in 
this estimate. There are examples where this estimate does not hold if the weights are not 
in ascending order. 

(2) At first glance one might claim that even the estimate 

™ 1^1 +n 

n 

n—1 

is valid. This is true in a lot of cases which we investigated. Nonetheless, there are 
examples where this estimate fails. 

The novel estimate 

™ -2- + n 

#Aw(g,TO) < 

n—1 


(SG Formula) 


n 
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is much sharper than the well established and widely used formula by Beged-Dov, cf. [3], 
(BD Formula) ^X^{q,m) < ^ 

n—1 

or the anisotropic tensor product estimate 


nw„ 


(TP Formula) 


#X^{q,m) < 


n—1 


+1 


Indeed, the new estimate implies (IBP Formula|) . This can easily be shown by induction: For 
m = 1, both estimates coincide. The induction step m i—>■ m + 1 follows with q = q + 
according to 


m+1 


m+1 


Y[i<l+ l|w||i) = {q + w„+i)™+^ = 


n—0 


“b 1 ) -n m+l-n 

q "'m+1 


> + (m + l)rwm+i 

> {q+{m + l)wm+i)q"' 


m+1 


> {q+{m + l)wm+i) n = nf? 


nw„ 


A numerical comparison of the exact number of indices, of our new bound and of the formula of 
Beged-Dov can be found in Figure [5] in the numerical examples. 

In practical applications, we will usually have to choose the level q. Hence, it is also interesting 
to examine how the computational cost behave, under the decay Assumption 14.21 with respect to 
the dimension m. 


Lemma 5.4. Let Assumption \4.2\ hold for some c > 1 and let m > 3. Then, we obtain that the 
number of indices in the anisotropic sparse grid on level q is bounded by 


(36) 


#Aw(g, m) < c(r) exp ( - log(log(TO)) ) = c(r) log(m)'^/ 


with a constant c(r) which is independent of m. 
Proof. From Lemma l5.21 we know that 


#Aw(g,m) < 


Next, we split the product into 


(37) 


n 


I = 


Wi 


-bl 


n—1 


q 

2 w2 


nw„ 




n—4. 


+ 1 • 


We estimate the last term by 

( 9 


n 


\nW'n 


+ i)=“p(Ei«g( + + i))s»p(E 


n=4 \ / ^n=4 

Due to Wn > log(n^), the sum in this estimate can be bounded by the following integral: 


E 


q 


< 


Js a;log(a:’') 


dx = 


r Js X log(x) 


dx 


^log(m) 2 ^ 


- dz = y log(log(m)) - log(log(3))). 


log(3) 
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The first three factors in dsa define a cubic polynomial in q and can thus be estimated by the 
exponential function according to 





< c(r) exp 


/log(log(3)) 


Hence, putting all together, we end up with 



< c(r)exp exp ^^(log(log(m)) - log(log(3)))^ 

< c(r) exp f - (log(log(m)) j . 


□ 


5.3. Convergence in terms of the number of quadrature points. The findings from the 
previous two sections can be summarized to an error estimate in terms of the sparse grid quadrature 
level q, that is 

(38) \{l-A^{q,m))U\ < C'i(5)e-«(i-^)||/^||c(s„), 

an error estimate in terms of the number of indices in the anisotropic sparse grid, 

(39) \{I- A^{q,m))fm\ < C2(w,/3)#Xw(g, 
an estimate of the computational cost, 

(40) N{q) := cost(.4w(g,ni)) < #X^{q,mf‘, 

and finally a sharp estimate on the number of indices in an anisotropic sparse grid, 


#X^{q,m) < 


f 9 


n—1 


\nWn 


1 


With these estimates at hand, we are now able to conduct the main result of this section. From 
(1391) and (I40p . we immediately obtain that the error in terms of number of quadrature points is 
bounded independently of the dimension m. 


Theorem 5.5. The error of the anisotropic sparse grid quadrature with Wn = log(Kn) can be 
bounded in terms of the number of quadrature points N(q) according to 

\{l- A^(q,m))fm\ < c(w,/3)iV(q')“(^“^^/^ 

for all fd < r. 


6. Numerical results 


This section is dedicated to numerical results in order to illustrate the theoretical findings. 
We will consider three different examples: At first, we consider a high dimensional quadrature 
problem, afterwards we have a look at a parametric diffusion problem as they usually occur in the 
context of uncertainty quantification and finally, we consider the approximation of quantities of 
interest from a diffusion problem on a random domain^ 


6.1. Pure quadrature problem. As a first example, we consider the quadrature problem 


/. p. 


/(y) := (o.6 + 0.2^n ’’y. 


for r = 2, 3,4. 


The dimension of the parameter is set to m = 1000. A reference solution is obtained by the 
anisotropic sparse grid quadrature on a higher level, featuring about 10® quadrature points. In 
order to validate these reference solutions, we have tested them against a quasi-Monte Carlo 
quadrature based on the Halton sequence. 


^The implementation of the sparse grid quadrature is available on https://github.com/muchip/SPQR. 
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Figure 3. Convergence of the quasi-Monte Carlo quadrature towards the refer¬ 
ence solution for the quadrature problem. 


As can be seen from Figure[31 the quasi-Monte Carlo quadrature converges towards our reference 
solution with a nearly linear rate for r = 3,4 and a slightly worse rate for r = 2. This behaviour 
is as expected from the results derived in |30) . 

Next, we consider the convergence of the anisotropic sparse grid quadrature. Obviously, in 
the absence of any decay, a genuine 1000-dimensional problem would be computationally not 
feasible. Thus, in order to determine the inherent dimensionality for each choice of the parameter 
r, we approximate the reference solution, which has been computed with m = 1000, by m = 
10 ,100,1000. 




Figure 4. Convergence of the anisotropic sparse grid quadrature with respect to 
the different parameter dimensions for r = 2 (left). Number of quadrature points, 
cardinality of the sparse index set and their estimates for dimension m = 100 
(right). 

The plot on the left-hand side of Figured shows the convergence for r = 2. It turns out, that we 
are able to recover the reference solution up to an error of order 10“^ by just considering to = 10 
dimensions. For to = 100, we already achieve convergence up to an error of order 10“® and, finally, 
for TO = 1000, we have convergence up to an error of order 10“^°. The observed convergence rate 
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is in all cases superlinear. Hence, the observed rate is much better than expected by Theorem 1 5.5 1 
since (/3 — l)/2 would be at most 1/2. 

The plot on the right-hand side of Figure |4] demonstrates that the bound on the number of 
quadrature points N(q) < is not sharp here. In fact, it seems that the number of 

quadrature points behaves more linearly in terms of the cardinality of the sparse grid index set. 
Indeed, for all the three different settings of this example (i.e., for r = 2,3,4), the number of 
quadrature points lies between the cardinality of the sparse index set and the novel upper bound 
(ISG Formulal) . This is the reason why we observe convergence rates that are better reflected by 
/3 — I than by (/3 — l)/2. 




Figure 5. Convergence of the anisotropic sparse grid quadrature with respect 
to different parameter dimensions for r = 3 (left). Number of points of the 
quadrature, cardinality of the sparse index set and their estimates for dimension 
m = 100 (right). 


Figure [5] depicts the situation for r = 3. Here for m = 10, we have already convergence up to 
an error of order 10“^. The choices m = 100,1000 both converge towards the reference solution 
up to an error of order 10“^^. For all choices of m, we observe an order of convergence that is 
greater than 2. The right-hand side of the figure shows that the number of quadrature points N(q) 
behaves very similar to log(m)'^/’' and that (ISG Formulal) is, overestimating N{q) by a factor 10 
for most of the considered values of q. 

Finally, on the left-hand side of Figure [6l we see the convergence of the sparse grid quadrature 
for r = 4. Here, m = 10 already provides convergence up to an error of order 10“®, whereas 
m = 100,1000 converge towards the reference solution up to an error of order 10“^^. We observe 
an order of convergence that is about 3 for all choices of m. On the right-hand side of Figure IHl 
we see that the number of quadrature points N(q) is bounded by 2 • ^X^{q, m) for the considered 
values of q. Moreover, it seems that the simplified bound from Lemma 15.41 with constant c(r) = 1, 
reflects the behaviour of #Xw((Z, w) quite well. 

6.2. Random diffusion problem. Let (H,J^, P) be a complete probability space and consider 
the diffusion equation 

(41) — dx{a{x,uj)dxu{x,u})) = 1 in = (0,1) for almost every a; £ H, 

where the random coefficient a(x, w) is uniformly bounded and elliptic meaning that 

0 < a := ess inf a(x,uj) and a := ess sup a{x,uj) < oo. 

{x,u:)GDxn {x,uj)eDxn 
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Figure 6 . Convergence of the anisotropic sparse grid quadrature with respect 
to different parameter dimensions for r = 4 (left). Number of points of the 
quadrature, cardinality of the sparse index set and their estimates for dimension 
771 = 100 (right). 


We arrive at a well-posed problem by complementing dni) with homogenous boundary conditions, 
i.e. 7 t( 0 ,a;) = u(l,a;) = 0 . 

The first step towards the solution for this class of problems is the parameterization of the 
stochastic parameter. To that end, one decomposes the diffusion coefficient with the aid of the 
Karhunen-Loeve expansion. Let the covariance kernel of a{x,oj) S Lp'[p,] {D)^ be defined by 

the positive semi-definite function 

C(x,x') '.= / (a(a:, w) — E[a](ai)) (a(a;', w) — E[a](ai')) dP(a;). 

Herein, the integral with respect to il has to be understood in terms of a Bochner integral, cf. HU. 
Now, let (Afc,(/ 7 fc) denote the eigenpairs obtained by solving the eigenproblem for the diffusion 
coefficient’s covariance, i.e. 

/ C{x,x')ipk{x')dx' = \k^k{x). 

Jo 

Then, the Karhunen-Loeve expansion of a(x, uj) is given by 

OO 

a(x,u;) = E[a](x) + ^ pX^<f„(x)Xn(uj), 

n—1 

where X„: ^ T C M for n = 1,2,... are centered, pairwise uncorrelated and L^-normalized 

random variables with Xn Note that the scaling factor -x/S stems from the Lp- 

normalization of the random varibles. We have additionally to assume that the random variables 
are independent. 

By substituting the random variables with their image, we arrive in the uniformly distributed 
case at the parameterized Karhunen-Loeve expansion 

OO 

a{x, Ip) = E[a](a::) + ^ \/K(pri{x)V3ipn, 

n—1 

where ipn G T. We define 7 „ = '/XP,\\^n\\L°°{D)- The decay of the sequence {'jn}n is important in 
order to determine the region of analytical extendability of the solution u, cf. Lemma l6.II 
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Truncating the respective Karhunen-Loeve expansion after m £ N terms, yields the parametric 
and truncated diffusion problem 

(42) - d^{am{x,y)do,Um{x,y)) = 1 in D = (0,1) for almost every y G T”". 

The impact of truncating the Karhunen-Loeve expansion on the solution is bounded by 

11'*^ ’^mW ~ C > 0, 

where £{m) —>■ 0 montonically as m —>■ oo, see e.g. miE]. Since the L^-norm is stronger than 
the L^-norm, this particularly implies the approximation estimate ([2]) for u and Um, where the 
modulus has to be replaced by the iLp (_D)-norm. 

Given the parametric solution Um{x,y), we are interested in determining properties of its dis¬ 
tribution. In our numerical examples, we focus on the computation of the solution’s moments. 
These are given by the Bochner integral 

:=^^Wm(a;,y)2“™dy. 

Especially, there holds E[um](a:) = Ml_^{x). 

We remark that the sparse grid quadrature straightforwardly extends to Bochner integrals, see 
e.g. m- This is due to the fact, that the Bochner integral corresponds for almost every x G (0,1) 
with the usual Lebesgue integral. 

It remains to provide the related regularity results that allow for an analytic extension of Um 
into the complex plane. The extendability is guaranteed by [3l Corollary 2.1], which has slightly 
been modified to fit our purposes. 


Lemma 6.1. The solution Um to ()42l) admits an analytic extension into the region for all r 
with 


I c —, where C(S) = k ^ ^ for arbitrary (5 > 0. 

C{S)k^+^'yk ^ 


Moreover, it holds that \\um\\c{i:m;H^(D)) bounded. 


Lemma [nm characterizes the region of analyticity and, therefore. Assumption 12 . ll for this appli¬ 
cation. The next lemma from [T] guarantees that Um '■ T"* —>■ Hq (D) also satisfies a one-dimensional 
error estimate for the polynomial approximation which gives rise to an estimate that is similar to 

(El- 


Lemma 6.2. Let X be a Banach space. Suppose that v G C{T\X) admits an analytic extension 
into El. Then, the error of the best approximation by polynomials of degree at most n can be 
bounded by 


(43) 
with Ki 


n + I+ t(. 


Thus, in view of ([4]), we end up with the estimate 

II(l(") - Q^”^)um(y*)||jji(^) < c(K„)exp( -log(K„)(2Af- l))||w„(y*)||c.(E„;ffl(_D)■ 

More generally we remark that by simply replacing all absolute values by the norms of the un¬ 
derlying Banach space X, the whole analysis we have presented so far for real valued functions 
directly transfers to Banach space valued functions fm'. T™ —>■ X which fulfill the estimate 

II (iW _ Q("))/™(y*)||^ < c(A„)exp(- log(K„)(27V- 1)) ||/m(y;^)||c(E„;A)- 


For this example, we employ two covariance kernels of the Matern class for :/ = 5/2 and v = 112, 
cf. [2T], i.e. 


1/ , y5p , 5p" 


C5/2(p) 7 ( 1 + — + ^ ) “ 
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Figure 7. Convergence of the anisotropic sparse grid quadrature for the random 
diffusion problem (u = 5/2 left, v = 7/2 right). 



Figure 8. Estimated number of points from the different formulas {v = h/2 left, 
V = 112 right). 


and 



Vlp 14p2 

~~r ^ 




where p := \x—x'\. The correlation length is in both cases set to £ = 1/2. The spatial discretization 
is performed with piecewise linear finite elements an a mesh with mesh size h = 2 “^'^, which 
results from 16384 equidistant sub-intervals. A numerical approximation to the Karhunen-Loeve 
expansion is computed by the pivoted Cholesky decomposition of the covariance operator with 
a trace error of e = 2“^®. This yields an approximation error of the underlying random field 
oi e = 2“^"*, see [15] for the details. The related truncation rank is given by m = 64 for C 5/2 
and m = 30 for C 7 / 2 . In addition, we set E[a](a;) = 2.5. From [TU Corollary 5], we know that 
7 „ < cn“® for C 5/2 and 7 „ < cn~'^ for C 7 / 2 . 

Since the solution of gu is not known analytically, we have again to provide a reference 
solution. The error with respect to the reference solution is measured relative to the norm of the 
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reference solution. This reference solution is computed by the quasi-Monte Carlo quadrature with 
Halton points and N = Ri 10® samples. 

For the anisotropic sparse grid quadrature, we choose the weights Wn according to Wn = log(K„) 
with Tn = l/7rf This would be the correct setting for a corresponding anisotropic tensor product 
quadrature. Hence, our anisotropic sparse grid quadrature is essentially a sparsification of the 
anisotropic tensor product quadrature, cf. m for more details on the anisotropic tensor product 
quadrature. To choose the same quantity r„ for the region of analyticity as for the tensor product 
quadrature seems to be a violation of Lemma l6.II Indeed, the assertion of this lemma is that the 
quantities t„, which describe the region of analytic extendability in each direction I]„, should be 
rescaled to = r„/(c(5)n^+^) in order to ensure analytic extendability into the tensor product 
domain S(t). Nonetheless, our experience suggests that the sparsification of the anisotropic tensor 
product quadrature yields an error which is nearly as good as the error of the anisotropic tensor 
product quadrature itself. 

For both cases v = 5/2 and v = 7/2, it turns out that we obtain similar convergence rates for 
all moments up to order four. The smoothness of the underlying covariance kernel has only little 
influence on the rate of convergence here. This might be caused by the fact that the eigenvalues 
in the Karhunen-Loeve expansion do not strictly decay in contrast to the previous example, but 
have some offset before the asymptotical rate is achieved. This is due to the small correlation 
length £ = 1/2. 

In addition to the convergence studies for the sparse anisotropic quadrature, we also provide 
results on the estimated number of indices contained in the sparse index set. To that end, we 
compare the tensor product estimate (ITP Formulal) and the formula by Beged-Dov (IBP Formula^ 
with the novel estimate proposed in this article (ISG Formulal) . It turns that the new estimate fits 
the cardinality of the sparse grid index set very well in comparison to the other estimates, cp. 
Figure El 



Figure 9. Different realizations of the random domain. 


6.3. Laplace equation on a random domain. As our third example, we consider the Laplacian 
on a random domain: 



The nominal domain Z?ref •= E[P] is given by an ellipse with semi-axis 0.3 and 0.2. The boundary 
perturbation is defined by means of a random vector field with V(dPref, w) = dDioj). This vector 
field is represented via a Karhunen-Loeve expansion according to 



where (/x S [0,27r) is the angle associated to x S 5Pref and A„ ^ W(—v^, v^). A visualization of 
several realizations of the random domain can be found in Figure E] By extending the boundary 
perturbation into space by a smooth blending function, it can be shown that the solution u to 
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the diffusion problem on the random domain exhibits a regularity similar to Lemma 16.11 with 
In '■= V%i||¥’nllvFLoo(£)_,^j.R 2 ), where the i~p^ correspond to the spatial functions of V’s Karhunen- 
Loeve expansion, see m for the details. 

As quantity of interest, we consider, likewise to the previous example, the first four moments of 
the solution, but this time restricted to the disc {x G : ||x ||2 < 0.05}. This disc is illustrated as 
grey shaded area in Figure |9l The numerical solution of the underlying boundary value problem 
is performed by an exponentially convergent boundary integral method based on collocation as 
proposed in [20) . The boundary is discretized by 200 points and a reference is again obtained by 

230 

quasi-Monte Carlo samples based on the Halton sequence. 



Figure 10. Convergence of the anisotropic sparse grid quadrature for the random domain. 

In Figure (TUI the relative errors of the moments computed by the sparse grid quadrature are 
depicted. Also in this application, we observe a rate of convergence which is greater than 1 for all 
moments under consideration and, therefore, again better than expected. 

7. Conclusion 

In the present article, we have considered the anisotropic sparse grid quadrature applied to a 
class of analytic functions. Under the assumption that the regions of analyticity for the particular 
dimensions are increasing algebraically, i.e. > cn’’ for some c, r G R+, we have derived dimension 
independent convergence of the anisotropic sparse grid quadrature. In order to estimate the cost 
of the quadrature, a novel estimate for the cardinality of the anisotropic sparse grid index set has 
been proven. Our numerical comparisons suggest that this novel estimate is much sharper than 
the widely used estimate by Beged-Dov. In particular, it has been shown that the estimate of 
Beged-Dov can easily be deduced from our novel estimate. Besides pure quadrature problems, the 
anisotropic sparse grid quadrature has been successfully applied to Bochner integrals which stem 
from the computation of the moments of elliptic diffusion problems with random coefficients or 
on random domains, respectively. 
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